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Surface states in topological insulators can be understood based on the well-known Shockley 
model, a one-dimensional tight-binding model with two atoms per elementary cell, connected via 
alternating tunneling amplitudes. We generalize the one-dimensional model to the three-dimensional 
case representing a sequence of layers connected via tunneling amplitudes t, which depend on the in- 
plane momentum p = {px,Py)- The Hamiltonian of the model is a 2 x 2 matrix with the off-diagonal 
element t{k,p) depending also on the out-of-plane momentum k. We show that the existence of the 
surface states depends on the complex function t{k,p). The surface states exist for those in-plane 
momenta p where the winding number of the function t(k,p) is non-zero when k is changed from 
to 27r. The sign of the winding number determines the sublattice on which the surface states 
are localized. The equation t{k,p) — defines a vortex line in the three-dimensional momentum 
space. Projection of the vortex line onto the space of the two-dimensional momentum p encircles the 
domain where the surface states exist. We illustrate how this approach works for a well-known model 
of a topological insulator on the diamond lattice. We find that different configurations of the vortex 
lines are responsible for the "weak" and "strong" topological insulator phases. The phase transition 
occurs when the vortex lines reconnect from spiral to circular form. We apply the Shockley model 
to Bi2Se3 and discuss applicability of a continuous approximation for the description of the surface 
states. We conclude that the tight-binding model gives a better description of the surface states. 

PACS numbers: 03.65.Vf, 73.20.-r, 31.15.aq 



I. INTRODUCTION 

Recent theoretical discovery plUH of topological in- 
sulators has stimulated active research in the field 
[l4*]. The key idea is that a time-reversal-invariant band 
Hamiltonian with a finite gap over the Brillouin zone 
(BZ) can be characterized by the topological indices Z2 
[2|-[j] . The topological indices distinguish trivial and non- 
trivial phases, which usually arise due to strong spin-orbit 
coupling in the system [l|, l6| . The topological Z2 indices 
are robust to moderate perturbations of the Hamiltonian 
and can change only if the energy gap is closed [ll-Q- 
Since the topological indices of vacuum and TI are dif- 
ferent ll2Ul5il. the boundary of TI should carry gapless 
modes [1^, [13 , similar to the chiral edge modes in the 
quantum Hall phases [I8] . Because of the bulk-boundary 
correspondence, the gapless surface states are topologi- 
cally protected against moderate perturbations. 

Theory of topological surface states has been studied 
in a number of works both in the tight-binding [19, 20] 
and continuous models j2llj26j. Many papers focused 
on the bulk-boundary correspondence, i.e., on proving 
that a sample with non-trivial topological numbers in 
the bulk should possess gapless excitations on the sur- 
face. The method of topological invariants, although be- 
ing very powerful, is often not physically transparent and 
not intuitive about the exact mechanism by which the 
topological numbers are related to the surface states. 

The purpose of our paper is to show that formation 
of the surface states in TIs can be understood based on 
the simple and well-known Shockley model [23, of 
the edge states. The Shockley model was also app lied 
to surface states in topological superconductors [29|, [30| ; 



however, we focus only on surfaces states in semicon- 
ductors. In Sec. [ni we review the one-dimensional (ID) 
Shockley model consisting of a chain of atoms connected 
via alternating tight-binding hopping amplitudes ti and 
^2- When a boundary is introduced in the system, e.g. 
by breaking the bond ^2, existence of the edge states is 
governed by the Shockley criterion. The edge state ex- 
ists if the greater tight-binding amplitude is broken at 
the boundary, i.e. if 1^2 1 > and it is localized on one 
sublattice. We show that the Shockley criterion can be 
formulated in terms of a topological winding number for 
the off-diagonal matrix element of the bulk Hamiltonian, 
thus connecting bulk properties with the surface states. 
In Sec, nil At we generalize the model to three dimensions 
(3D) by replacing atoms by the two-dimensional (2D) lay- 
ers parallel to the xy plane and assigning the in-plane 
momentum dependence p = {px^Py) to the interlayer 
tunneling amplitudes ti and ^2- In Sec. IIIIBl we study 
vortex lines in the 3D momentum space [31-33], where 
the off-diagonal matrix element of the bulk Hamiltonian 
vanishes. We show that the projection of the 3D vortex 
lines onto the 2D in-plane momentum space encircles the 
domain where the surface states exist. We observe that 
the tight-binding TI Hamiltonians [1, H, i, IH have the 
Shockley-model structure and can be understood using 
our approach. In Sec. IIV At we illustrate the Shockley 
mechanism for the Fu-Kane-Mele model on the diamond 
lattice [sf. We show how the surface states evolve when 
the parameters of the Hamiltonian vary. In Sec. IIVB( 
we show that reconnection of the vortex lines represents 
a phase transition in the TI Hamiltonian. The spiral vor- 
tex lines correspond to a phase with an even number of 
Dirac cones (the "weak" TI phase), while the circular 
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FIG. 1: (Color online) Panel (a): ID chain of atoms with alternating tunneling amplitudes ti and t2 representing the Shockley 
model, Eq. Panel (b): The bulk energy spectrum of the system, Eq. (|8]), with a non-zero gap for |ti| ^ \t2\. Panel (c): The 
exponentially decaying edge state, Eq. ()18p . for |ti|/|t2| < 1 with the penetration depth ^ = l/ln|t2/ti|. 



vortex lines correspond to a phase with an odd number 
of Dirac cones (the "strong" TI phase). In Sec. |Vl we 
apply the Shockley model to describe the surface states 
in Bi2Se3, which is formed by the quintuple layers of Bi 
and Se j2l|, HI, The electronic structure of this 

material near the Fermi level can be well described by 
the hybridized Pz orbitals located near the outer layers 
of the quintuplets [H, [sj- Thus, the Shockley model 
with the intra- quintuplet and inter-quintuplet tunneling 
amplitudes ti ant t2 gives a plausible description of this 
material. Surface states have complementary properties 
depending on how the crystal is terminated [3^. Break- 
ing the t2 amplitude introduces a cut between the quin- 
tuplets. In this case, the surface states have a Dirac 
cone in the Brillouin zone (BZ) center [35|, l36|. Breaking 
ti introduces a cut inside the quintuplet. In this case, 
the Shockley model predicts the surface states with the 
Dirac cones on the boundary of the BZ. The similar ef- 
fect was considered for the Bii-^^Sb^; alloy in Ref. [38]. 
In Sec. lVB) we discuss whether a continuous approxima- 
tion for the TI Hamiltonian gives a good description of 
the surface states. We conclude that the tight-binding 
models are better suitable for the description of the sur- 
face states. Then, in Sec. IVH we generalize the Shockley 
model by including additional tight-binding amplitudes. 
For all these models, we find that the edge state is always 
localized on one sublattice, which is rarely mentioned in 
the TI literature. 



II. ID SHOCKLEY MODEL 

A. The original Shockley model 

In this section, we briefly review the Shockley model 
[27I , [28| and its properties. Let us consider a ID linear 
chain of atoms shown in FiglTJa). The unit cell contains 
two atoms labeled as A and B, which are connected via 
the alternating nearest-neighbor complex tight-binding 



amplitudes ti and t2- So, the Hamiltonian of the model 
is 

H = Y. ^t(^) [u^j{z) + V^[z - 1) + V^^[z + 1)] ,(1) 

Here, z is the integer coordinate of the unit cell, ti and ^2 
are the intra-cell and the inter-cell tunneling amplitudes, 
and ^{z) is the spinor 



(3) 



where tl^ai^) and tl^bi^) are the wave functions on the 
sites A and B. In the Fourier representation ^(^) = 
/q g e'^^ ^(fc), the Hamiltonian is 

H= ^^If^{k)H{k)^{k), (4) 

Jo 



where 

H{k) = U + Ve-^'^ + yte"= = ^ « ^ ^ (5) 

is a 2 X 2 matrix acting in the AB sublattice space, and 
t{k) = ti + t2e'^ = + t2q, q = e'\ (6) 
Then, the Schrodinger equation 



e{k)\f^a 

m i I 



E 



(7) 



gives two particle-hole symmetric energy bands with the 
eigenvalues E{k) and eigenfunctions ^(/c) 



E{k) = ±m\, 

giarg[^(fe)] 



±1 



(8) 
(9) 



The energy spectrum has a gap if |ti | ^ |t2 1, as illustrated 
in Fig. for real ti and ^2- Notice that the bulk wave 
function ([9]) has equal probabilities on both sublattices. 
In contrast, as we shall see below, the wave function of 
an edge state is localized only on one sublattice. 

A boundary to the ID lattice can be introduced by 
cutting either ti or t2 link. Let us consider a half-infinite 
system forz>l, z = l,2,3, corresponding to the 
cut of the t2 link. In this case, the atom A is exposed 
on the edge, as shown in Fig. [U^a). Mathematically, the 
boundary condition is introduced by requiring that the 
wave function vanishes at the fictitious site z = and at 
infinity 

^a(o) = o, Mo) = o, (10) 

^a(+Oo)=0, M + ^) = 0. (11) 

It is shown in Appendix El that the edge state can exist 
only for E = 0. So, we substitute E = into Eq. ([7]) and 
find that the wave functions on the A and B sublattices 
decouple 



t{k) -- 



ik \ 



t2e 
Vtle 



0, 



0, 



(12) 
(13) 



where k is now a complex wave- number, so t*(k) is not 
a complex conjugate of t{k). In the real space, Eq. (p!3|) 
can be written as a recursion relation 







(14) 



for z > 1. Using this recursion relation and the boundary 
condition ^^^(O) = 0, we find that il)h{z) vanishes for z > 
1. In contrast, the real-space representation of Eq. ([12]) 



llja{z)h ^ t/jaiz ^ l)t2 







(15) 



for 2: > 1 does not involve '0a(O) from Eq. (p^Oj) . So, the 
solution on the A sublattice is 



(16) 



where, qo is obtained by solving the equation t{ko) = 0, 
following from Eq. ([12]) for a complex wave-number ko 



Qo 



(17) 



Depending on whether |go| < 1 or \qo\ > 1, the solution 
in Eq. ([T6|) either satisfies the condition ([TT]) at infinity 
or not. If 1^2! > 1^1 1, then \qo\ < 1, as shown in Fig. |2]^a), 
and the wave function ([T6]) exponentially decays at z ^ 
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FIG. 2: (Color online) Topological formulation of the Shock- 
ley criterion (|19p . Panels (a) and (b) compare the two cases, 
where the root ^0 (the red dot) lies inside or outside the unit 
k G (0,27r)}. An edge state exists for 
and does not exist for |^o| > 1, panel 
(b) . An alternative formulation in terms of the winding num- 
ber ()20|) is illustrated in panels (c) and (d). The edge state 
exists if the winding number is non-zero, panel (c), and does 
not exist if the winding number is zero, panel (d). 



+00, as shown in Fig. Hj^c), so the edge state exists. In 
contrast, if I ti I > |t2|, then |go| > 1, as shown in Fig.|2]^b), 
and the wave function ([T6|) exponentially grows at z ^ 
+00, so an edge state does not exist. To summarize, by 
solving Eqs. ([T2]) and ([T3]) with the appropriate boundary 
conditions ([TQ|) and ([TT]) , we obtain the zero-energy edge 
state 



^o(^) 



which exists only if 



1901 



9o" 



Eo = 0. 



\t2\ 



(18) 



(19) 



Equation ([T9|) constitutes the Shockley Criterion: In the 
ID tight-binding model with alternating tunneling ampli- 
tudes given by Hamiltonian (Q])^ the edge state exists if 
the bond of the greater magnitude is broken at the bound- 
ary. 

Let us now discuss an alternative formulation of the 
Shockley criterion ([T9|) in terms of the winding number 
W defined as 



W = 



1 

2^i 



dk-^\nt{k). 
dk 



(20) 



The winding number W represents the phase change 
of the complex function t{k) when the real variable k 



changes from to 27r. The function t{k) also defines a 
closed contour 



C' = {t{k) =ti+ ^2e*^ k e (0, 2tt)} 



(21) 



in the 2D plane of (Re t,Imt), as shown in Fig. [21 panels 
(c) and (d). If |(7o| < 1, or equivalently \t2\ > |ti|, the 
contour winds around the origin (red dot), as shown 
in panel (c). If \qo\ > or equivalently \ti\ > 1^2 1, the 
contour does not wind around the origin, as shown in 
panel (d). So, the Shockley criterion (p!9|) can be formu- 
lated in terms of the winding number 



W 



1, edge state exists, 

0, edge state does not exist. 



(22) 



This formulation will be useful for the discussion of the 
generalized 3D Shockley-like models in the following sec- 
tions. 
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FIG. 3: (Color online) 3D generalization of the Shockley 
model described by Hamiltonian ()28p with h(p) defined by 
Eq. (|3T]) . The arrows show the staggered direction of the 
Rashba vector n. 



B. On-site energies in the ID Shockley model 



III. 3D SHOCKLEY-LIKE MODEL 



Let us further generalize the model and include on-site 
energies Sa and £5 in Hamiltonian ^ 



A. Generalization to the 3D case 



H{k) 



( Sa t*{k) 
\ t{k) Sb 



(23) 



As shown in Eq. (p!8|) for £^=£5 = 0, the edge state so- 
lution is localized on the A sublattice. Therefore, adding 
the on-site energy Sa simply shifts the energy of the edge 
state without changing its wave function irrespective of 
£i). So, if criterion (fT9|) is satisfied, the edge state is lo- 
calized on the A sublattice and has the energy 



^0 = ^a- 



(24) 



It is also convenient to transform the Hamiltonian to the 
symmetrized form 



Hik) = 



The offset (£:a+£:5)/2 just uniformly shifts all energies and 
will be omitted in the rest of the paper, so the Hamilto- 
nian becomes 



H{k) 



h t*{k) 
t{k) -h 



(26) 



The bulk spectrum of the Hamiltonian (|26]) is generally 
gapped 



Let us generalize Hamiltonian (j26|) to the 3D case. In- 
stead of alternating atomic sites, let us consider a se- 
quence of alternating layers A and B perpendicular to 
the z direction, as shown on Fig. [3l Now, all parameters 
of Hamiltonian (|26]) acquire dependence on the in-plane 
momentum p = (px^Py) 



H{k) = 



Hp) t*{k,p) 
t{k,p) -Hp) 



The off-diagonal matrix element 

t{k,p)=h{p)+t2{p)e 



ik 



(28) 



(29) 



describes the p-dependent inter-layer tunneling ampli- 
tudes, while h{p) represents the intra-layer Hamiltonian. 
Throughout this paper, we denote the in-plane momen- 
tum as p = {px-tPy) and the out-of-plane momentum in 
the z direction as k [39]. 

For a fixed value of the in-plane momentum p, Hamil- 
tonian ([28|) reduces to the ID model (j26j), for which the 
edge state were studied in Sec. [Ill The surface states ex- 
ist for those in-plane momenta p where criterion (p!9|) is 
satisfied. The surface states are localized on the A sub- 
lattice, and the energy spectrum Eo{p) of the surface 
states is determined by the in-plane Hamiltonian h{p) 



Eo{p) = Hp). 



(30) 



E{k) = ±^/h^^\t{kj^. (27) 



Now let us illustrate how to obtain a linear Dirac-like 
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FIG. 4: (Color online) Energy spectrum of the 3D Shockley 
model described by Hamiltonian (|28p in the vicinity of p = 0. 
The spectrum of the bulk states, Eq. (|33p , is shown by the 
solid parabolas in both panels. According to the Shockley 
criterion, surface states exist if < |t2|, panel (a), and do 
not exist otherwise, panel (b). The surface states have the 
linear dispersion, Eq. (|32|) . shown by the transparent Dirac 
cone in panel (b). 




FIG. 5: (Color online) The thick blue curve is a vortex line in 
the 3D momentum space defined by Eq. (|34p . Its projection 
onto the 2D momentum space p defines the boundary of the 
shaded area, where the surface states exist. 



B. Vortex lines in 3D momentum space 



energy dispersion for the surface states in model (|28|) . 
Let us consider h{p) in the vicinity of a time-reversal- 
invariant momentum point p = Pi, e.g. near p = 0. In 
the presence of the spin-orbit coupling, a generic Hamil- 
tonian h{p) must be bilinear in the spin Pauli matrices 
a and momentum p to satisfy the time-reversal invari- 
ance. For example, h{p) can have the Rashba spin-orbit 
coupling form 

h{p) = v{axPy - (JyPx) = v{a x p) ■ z, (31) 

where v has the dimension of velocity. So, h{p) is now 
a 2 X 2 matrix and H{k) a 4 x 4 matrix. Notice that 
the diagonal term ^h{p) in Hamiltonian (f28|) has oppo- 
site signs on the A and B sublattices. This corresponds 
to staggered direction of the Rashba vector n = ±i on 
different layers for the spin-orbit coupling vn{a x p) as 
shown in Fig. [3l The surface states exist only if < 1^2 1, 
and the spectrum of the surface states has linear depen- 
dence on |p| 

Eo{p) = ±v\pl (32) 

as shown in panel (a) of Fig.lH The wave functions of the 
surface states have in-plane spin-polarization perpendic- 
ular to the momentum p. On the other hand, the bulk 
spectrum is parabolic in the vicinity of p = 

E\k,p) = \t{k,0)\^ + v^p\ (33) 

as shown in both panels of Fig. |4]by solid colors. 



In principle, the tunneling amplitudes ti{p) and t2{p) 
may depend on the in-plane momentum p. So, the sur- 
face state existence criterion can only be satisfied in a 
certain domain of the 2D momentum space p. In this 
section, we discuss how to identify this domain for the 
Hamiltonian (|28|) . 

Let us consider the equation 

t{k,p) = (34) 

for the complex function t{k^p) in Eq. (|29|) . It is equiva- 
lent to two equations Re t{k,p) = and Im t{k^p) = 0, 
which define a line in the 3D momentum space (k^p). In 
general, the complex- valued function t{k,p) has a phase 
circulation around the line where it vanishes, i.e. Eq. (|34|) 
defines a vortex line in the 3D momentum space [3L-.33J . 
As an example, such a vortex line and its projection on 
the 2D momentum space p are shown in Fig. O Phase 
winding of the function t{k^p) along an arbitrary contour 
7 can be calculated as 

where the notation / = (/c,p) is used for brevity. For 
instance, the phase winding along the contour 73 around 
the vortex line in Fig. [5] is non-zero 

W{j3) = 1. (36) 

Because the BZ is periodic in /c, we can also define a 
closed contour by varying < /c < 27r for a fixed value 
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FIG. 6: (Color online) Illustration of the diamond crystal 
structure and the tight-binding model described by Hamilto- 
nian (|39p . The lattice has two atoms in a unit cell shown by 
the red (A) and blue (B) spheres. 



of the in-plane momentum p. Such contours 71 and 72 
are shown in Fig. [5l and the phase windings (|35]) are well 
defined for these contours. The contours 71 and —72 can 
be merged into the contour 73. So, the following equation 
holds 

1^(73) = W(7i)- ^(72). (37) 

Given Eq. (j36l) and the condition (j22j) that 1^(71,2) > 0, 
we find that the winding numbers are W{ji) = 1 and 
^(72) =0. Since a non-zero winding number is required 
for existence of the surface states according to Eq. (|22]) . 
we conclude that the surface states exist for the 2D mo- 
menta p in the shaded area of Fig. [5] and do not exist 
outside. Thus, we have shown that the projection of the 
vortex line (|34|) onto the 2D momentum space p defines 
the domain where the surface states exist. 



IV. DIAMOND MODEL 
A. Hamiltonian and surface states 

In this section, we illustrate how the Shockley model 
can be applied to study the surface states for a particular 
TI model of Ref. [S]. However similar approach can be 
applied to other models [l, [HI, [TqI . 

Let us consider a tight-binding model on the diamond 
lattice shown in Fig. [6l The diamond lattice has two 
equivalent atom positions denoted by A (red) and B 
(blue). Atoms of each type form 2D triangular lattices, 
so that the A and B layers alternate along the z direc- 
tion similarly to Fig. [3l The nearest A and B layers 
form a distorted graphene lattice. So, when viewed along 




FIG. 7: (Color online) Lines of constant value for the 
graphene spectrum function |ti(p)|/|ti| = C, for C — 
0.5, 1, 2, within the Brillouin zone (BZ), denoted by the 
dashed lines. The contour lines degenerate to points at the BZ 
corners (K and points) at C = and at the BZ center at 
C = 3. The thick red dots denote the time- reversal- invariant 
momenta points (|45|). 



the z direction, the structure looks like the ABC-stacked 
graphite lattice. We define the nearest-neighbor vectors 
an, n = 1, 2, 3, 4, as shown in Fig.[6l as well as the vectors 

^i=a3-a2 = (l/2,-V^/2), 
^2 = ai -as = (1/2,^3/2), (38) 
^3 = 02 - ai = (-1 , 0), 

which are the in-plane elementary translation vectors of 
the unit length \5n\ = 1. 

The Hamiltonian of the model has the form of Eq. (|28|) 

where the unit cell consists of the A and B atoms con- 
nected by the vector as. The off-diagonal part 

t{k,p)=h{p)^t2e'\ (40) 
ti(p) =ti(l + e-*P^^ +e*^^^), (41) 

describes the nearest-neighbor tunneling between the A 
and B sublattices with the amplitude ti along the vectors 
Ojj, n = 1,2,3, and the amplitude ^2 along the vector 
[4^. In Eqs. (go]) and (gT]) . we distinguish between the 
in-plane- momentum-dependent function ti{p) and the 
tight-binding amplitude ti. Equation (|^T]) describes the 
well-known tight-binding spectrum of graphene [41^ 

\ti{p)\/\ti\ = V3 + 2 cos(p^i) + 2 cos(p^2) + 2 cos(p^3) 

(42) 
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FIG. 8: (Color online) The plot of the particle-hole symmetric 
spectrum Eo(p), Eq. (|44p . induced by the spin-orbit Hamilto- 
nian h{p) ()43p . In the vicinity of the time reversal invariant 
points, shown with the thick red dots, the Hamiltonian ()43|) 
becomes linear in momentum. The dashed line denotes the 
boundary of the BZ. 



The contour plots of \ti{p)\/\ti\ = C, for C = 0.5, 1, 2, 
are shown in Fig. [71 Note that ti (p) has the linear Dirac- 
like dependence on the momentum p at the BZ corners, 
K and points in Fig. [71 



The diagonal term h{p) in Hamiltonian 
the spin-orbit interaction ^] 



describes 



Hp) 



2V2 



Aso ^ ^iji {o- • [ai X aj]) sm{pSi), 

i,j,^=l,2,3 



(43) 

where Aso is the strength of the spin-orbit coupling, and 
eiji is the antisymmetric tensor. For simplicity, we do not 
include the inter-layer spin-orbit coupling involving the 
vector a4 in Hamiltonian (|43|), unlike in Ref. [3]. Hamil- 
tonian (|43l) has a gapless particle-hole symmetric spec- 
trum Eo{p) = ±^W(p), 

^o(p)/Aio = h\p)/Klo = sin^(p^i) + (44) 
= sin^(p^2) + sin^(p53) + sin(p^i) sin(p^2) + 
+ sin(p^i) sin(p^3) + sin(p^2) sin(p^3), 

which is shown in Fig. [H The energy Eo{p) vanishes at 
the four time-reversal-invariant momenta (TRIM) 



G {r,Mi,M2,M3}, 



(45) 



where F is the BZ center, and Mi, M2, M3 are at the 
centers of the BZ boundary, as shown in Figs. [71 and [H 
Hamiltonian (1^3]) is bilinear in the momentum and spin 





t2 broken 


ti broken 


TI 


0< |t2| < \ti\ 




Mi,2,3, r 


Weak 


\tl\ < \t2\ <3\h\ 


Mi,2,3 


F 


Strong 


3|ti| < \t2\ 


Mi,2,3, r 




Weak 



TABLE I: The table shows the points in the BZ where the 
surface states exist depending on the parameters of the model 
and which bond is broken at the surface. According to Fig. [H 
the surface states have the Dirac cones at the corresponding 
points. Letters Mi, M2, M3, F denote positions of the TRIM 
points (|45)) . 



operators in the vicinity of the TRIM points 



Aso 
Aso 



V3 1 



(46) 



2V2 

7f 



(47) 



Thus, the energy spectrum Eo{p) has the shape of the 
Dirac cones in the vicinity of TRIM points, as shown in 
Fig. [H It is important to distinguish the linear, Dirac- 
like, behavior of the off-diagonal term ti (p) in the vicinity 
of the BZ corners (K and points) and of the diagonal 
term h{p) in the vicinity of the TRIM points, which are 
different sets of points in the BZ. 



The bulk spectrum of Hamiltonian (|39l) 

E^{k,p) = \t{k,p)\^ + El{p), 



(48) 



contains contributions from both the diagonal h{p) and 
the off-diagonal t{k,p) terms. The bulk spectrum be- 
comes gapless when both contributions vanish for some 
momenta {k,p) 



Eoip) ■- 
t{k,p) 



0, 
-- 0. 



Given Eq. (HOJ, Eq. (jSOl) is equivalent to 

\tl{p)\ = \t2\, 



(49) 
(50) 



(51) 



which defines a contour line in Fig. [71 Condi- 
tions (|49)) and (|5Ql) can be satisfied simultaneously only 
for special values of the parameters ti and ^2- The 
bulk spectrum becomes gapless, for \t2\ = when 
the contour Hue ([5T]) passes through the TRIM points 
Ml, M2, M3, and for \t2 \ = 3|ti|, when it passes through 
the F point. 

Hamiltonian (|39l) has the Shockley form, Eq. (|28]) . 
Therefore all the conclusions of Sees. [Til and IIIII apply 
here, including the criterion (p^Qj) for existence of the sur- 
face states. We find that the surface states have the 
dispersion Eq (p) and exist for those in-plane momenta p 
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where the following condition is satisfied 

1*1 (P) I < l*2|. 



(52) 



The boundary of this domain is given by Eq. (pT]) . When 
we change the parameter t2 while keeping ti fixed, the 
Hamiltonian undergoes a transition between the phases 
with odd and even numbers of surface Dirac cones, called 
the "strong" and "weak" TI phases in Ref. [3] . For small 
1^2! ^ 1^1 1 5 the contour lines given by Eq. (|5T]) wind 
around the BZ corners K and (see Fig. [71 for C = 0.5), 
and criterion (|52]) is satisfied in the small area inside. The 
surface states do not include the Dirac cones of Eo{p)^ 
shown in Fig. [51 because the TRIM points (red dots) are 
in the area where Eq. (|52]) is not satisfied. Thus, Hamil- 
tonian (|39|) is in the "weak" TI phase in this case. For 
1^1 1 = 1^2!, the contours ([5T]) become straight lines pass- 
ing through the TRIM points Mi, M2, M3, as shown in 
Fig.[7]for C = 1. So, both Eqs. (HI]) and dSH are satisfied 
at the TRIM points, and the bulk spectrum, Eq. (|48|) . be- 
comes gapless. This marks a transition to the "strong" 
TI phase. When \ti\ < \t2\ < 3|ti|, the contour forms a 
circle around the BZ center, see Fig. [7]for C = 2. Cri- 
terion Eq. ([52]) is satisfied in the exterior of the circle, 
and so the surface states contain the Dirac cones at the 
TRIM points Mi, M2 and M3. When ^2 reaches the crit- 
ical value 1^2! = 3 1^1 1, the contour ([5T|) shrinks to the 
single point F. The bulk spectrum becomes gapless, and 
this marks a transition to the "weak" TI phase again. 
For 1^2! > 3 1^1 1, the Shockley criterion is satisfied every- 
where in the BZ, so the surface states include the Dirac 
cones for all TRIM points ([45]) . 

As discussed above, the Shockley criterion ([52]) is writ- 
ten for the case where the ^2 bond is broken at the surface. 
If, on the other hand, the crystal termination is such that 
the ti bond is broken at the surface, the existence crite- 
rion for the surface state becomes complementary to the 
criterion ([52l) 




1*1 (P) I > 1*21. 



(53) 



So, the surface states now exist for those momenta p 
where they did not exist in the case of the broken bond 
t2 and have the Dirac cones at the complementary TRIM 
points. This is summarized in Table [H which shows the 
Dirac cones belonging to the surface states depending on 
whether ti or ^2 is broken at the surface. In the "strong" 
TI phase, there is an odd number of the Dirac cones in the 
surface states, so, at least, one surface Dirac cone always 
exists. In contrast, in the "weak" TI phase, there is an 
even number of the surface Dirac cones, so the surface 
states may disappear under certain conditions. 



B. 3D vortex lines 



FIG. 9: (Color online) Vortex lines in the 3D momentum 
space defined by Eq. (|54p . and shown for different values 
of the parameters: (a) |ti| > |t2|, (b) |ti| = |t2|, and (c) 
\ti\ < \t2\ < 3|ti|. The vortex lines are shown by the thick 
lines with the arrows representing vorticity. The thin lines 
show projections of the vortex lines, which encircle the shaded 
area in the 2D momentum space p, where the Shockley crite- 
rion ([52]) is satisfied and the surface states exist. The dashed 
lines show the boundaries of the BZ. The part of the vortex 
lines residing in the first BZ is highlighted in red in panel 
(b). The three panels show the evolution of the vortex lines 
with the change of the parameters of the Hamiltonian. At 
|ti| = |t2|, the vortex lines reconnect at the TRIM points and 
change their topology from spirals for |ti| > \t2\ to the loop 
for |ti| < \t2\. The change of the vortex lines topology is 
responsible for a transition from the "weak" to "strong" TI 
phase in the Hamiltonian (|39p . 



In the previous section, we showed that the 2D con- 
tour defined by Eq. ([5T]) represents the boundary sepa- 
rating the domain in the 2D momentum space where the 



Shockley criterion is satisfied. On the other hand, the 
contour ([5T]) is just the projection of the 3D vortex Hue, 
defined by Eq. ([50]) 



(54) 



onto the 2D momentum space, as discussed in Sec. IIIII 
Let us discuss evolution of these 3D vortex hues with 
the change of the parameters ti and ^2- In the vicinity 
of the BZ corners po = (±47r/3,0), K and points in 
Fig. [71 where the function ti{p) vanishes, Eq. (|4T]) can 
be hnearized 



hiPo^p) 



V3 



(55) 



So, for 1^2 1 < l^il, Eq. §^ with ti{p) defined in Eq. (|55 
describes spirals in the 3D momentum space {k,p) ^] 




2 ^2 / ± cos k 
\^ti I sink 



(56) 



as shown in Fig. [9l^a). Projections of these spirals onto 
the 2D momentum space p encircle the corners K and 
of the 2D BZ. With the increase of ^2, the spirals 
grow until t2 reaches the critical value 1^2 1 = At this 
point, the vortex lines reconnect as shown in Fig. [9jb) 
and transform into three families of straight lines ob- 
tained by intersections of the planes 

{p6i = TT + 27rn} f]{p6s = -k ^ 27rm}, (57) 
{pS2 = TT + 27rn} f]{pS3 = /c + 27rm}, (58) 
{pSs = TT + 27rn} f^{k = tt + 27rm}, (59) 

where n and m are independent integers. The part of 
these lines residing in the first BZ forms a loop high- 
lighted in red for clarity in Fig.[9l^b). With the further in- 
crease of ^2, the vortex line detaches from the BZ bound- 
ary and becomes a closed loop, as shown in Fig.[9jc). In 
the vicinity of the F point, the function ti{p) given by 
Eq. (|4T]) can be expanded to the second order in p, so 
the vortex line defined by Eq. ([5^ is given by the inter- 
section of a cylinder and a plane in the 3D momentum 
space (A:,p) 



t2 



(60) 

For the critical value \t2\ = 3|ti|, the vortex line shrinks 
to the F point and then disappears for 1^2 1 > S\ti\. 

So, we observe that the vortex lines change their 
topology at the critical values of the model parameters 
1^2! = 1^1 1 and 1^2! = 3|ti|. These are the critical values 
where the transitions happen between the "weak" and 
"strong" TI phases. So, the configuration of the vortex 
lines (|54|) is directly related to the topological phase of 
the fuh Hamiltonian Eq. (|39|). 




FIG. 10: (Color online) The crystal of Bi2Se3 is formed by 
quintuple layers, schematically shown by the blue boxes. Each 
quintuplet consists of the alternating layers Se-Bi-Se-Bi-Se. 
The tight-binding tunneling amplitudes ti and t2 connect the 
orbit als of the outermost edges of the quintuplets. Then, de- 
pending on whether t2 or ti is broken at the surface, as shown 
by the red line in panels (a) and (b), surface states occur in 
different regions of the 2D momentum space, as shown in 
Panels (c) and (d). 



V. SHOCKLEY MODEL DESCRIPTION OF 

BiaSes 

A. General analysis 

Despite its simplicity, the Shockley model may be di- 
rectly relevant to the description of real materials, such 
as Bi2Se3. The crystal of Bi2Se3 is formed by a se- 
quence of quintuple layers [211, [H, [H, S, [s^- Each 
quintuplet consists of five alternating layers of Bi and 
Se, as sketched in Fig. [TOTa). Chemical bonding within 
the quintuplets is relatively strong, whereas the inter- 
quintuplet Van der Waals attraction is relatively weak. 
So, the natural cleavage plane lies between the quintu- 
plets, as shown in Fig. [TOTa). 

For the relevant energy interval near the Fermi level, 
the electronic structure can be captured by considering 
the electronic orbitals localized near the outermost layers 
of Se within the quintuplets [s^, as shown by the thick 
lines in Fig. [TOT a). Then, the Shockley amplitudes ti 
and t2 describe the intra- and inter-quintuplet tunneling 
between these orbitals, as shown in Fig. [TOTa). The tun- 
neling amplitudes ti and t2 may depend on the in-plane 
momentum p. 

As shown in the previous section, the Shockley surface 
states strongly depend on how the crystal is terminated. 
When the crystal is cut between the quintuplets, and t2 
is broken on the surface as shown Fig. [TOf a) and realized 
experimentally, a single Dirac cone is observed at the 
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BZ center [35|, as shown in Fig. [TOTc). So, in terms of 
the Shockley model, the surface state existence criterion 
1^1 (p) I < 1^2 (p) I is satisfied at the BZ center and not 
satisfied at the BZ boundary. 

In principle, the surface can also be introduced by cut- 
ting the quintuplet layer and breaking the bond ti, as 
shown in Fig. [ToT b). To the best of our knowledge, this 
type of surface has not been observed in Bi2Se3. In the 
previous section, we found that, for alternative crystal 
terminations, the surface states have the Dirac cones at 
the complementary TRIM points of the 2D BZ. Thus, 
we conclude that, when the quintuplet is broken at the 
surface as in Fig. [TOT b). the surface states should have 
Dirac cones at the boundary of the 2D BZ as shown 
in Fig. [TOT d). A similar prediction was made for the 
BiajSbi_a; alloy in Ref. [ssj. 

We can estimate the Shockley tunneling amplitudes 
ti{p) and t2{p) for p close to the F point based on 
the band-structure calculations of Ref. [21]. As dis- 
cussed in Sec. Ill A| the extreme values of the energy 
gap can be obtained from the off-diagonal matrix el- 
ement = ti + ^26*^ ^ = ti ± t2. We compare 
these values with the band structure along the direction 
k G (0,7r) for the fixed in-plane momentum p = 0, which 
is shown in Fig. 2(b) of Ref. [21]. From the set of equa- 
tions ti-\-t2 = 0.28 eV and h —t2 = —0.6 eV, we obtain 
the following estimate 

h = -0.16 eV, t2 = 0.44 eV. (61) 

Note that the geometric distance between the orbitals 
on the adjacent quintuplets is shorter than the distance 
between the orbitals within the quintuplet. Therefore, in 
the vicinity of the F point, the inter-quintuplet tunneling 
1^2! should be greater than the intra-quintuplet tunneling 
|ti|, which is consistent with Eq. ([6T]) . 



B. Continuous approximation 

In previous sections, we have shown that, in the Shock- 
ley model, existence of the surface states relies explicitly 
on the tight-binding nature of the model. However, con- 
tinuous models [21 --26] are also widely used to describe 
the surface states in the TI models and in real materials, 
such as Bi2Se3. A continuous approximation is obtained 
by expanding the Hamiltonian in the powers of the mo- 
mentum k in the z direction. This is equivalent to disre- 
garding the BZ periodicity for the momentum k and tak- 
ing the limit where the size of the elementary cell in the 
z direction goes to zero. In this subsection, we examine 
the applicability of the continuous-limit approximation. 



diagonal elements vanish 

t{k)=h^t2e'\ (63) 

Without loss of generality, let us make an assumption 
that t2 > 0. If the energy gap \t{k)\ reaches minimum at 
A: = 0, then ti < as in Eq. (|6T]) . Then, we expand t{k) 
to the first order in k around k = 

t{k) =ti^t2^it2k. (64) 

The tight-binding boundary conditions (p!Q|) and (pT]) cor- 
respond the following boundary conditions [34] for the 
continuous approximation 

V^a(^^0)/0, V^5(^^0)=0, (65) 
ilja,b{z ^ 00) = 0. (66) 

Using these boundary conditions, we solve the 
Schrodinger equation H"^ = for Hamiltonian ([62|) with 
the continuous t(/c), Eq. (|64|) . and obtain the surface state 

*o(^) = J j e*°^ (67) 

The exponential decay length in Eq. ([67|) is given by the 
parameter 

^0 = ^ (^1 + ^) , (68) 

which is the root of the equation t{ko) = ti -\- 12 -\- it2ko = 
0. Boundary conditions (|66|) are satisfied if Im /cq > or 
equivalently ti > — ^2; otherwise, the surface state does 
not exist if ti < — ^2- So, the continuous model (|64|) with 
the appropriate boundary conditions (|65]) and (|66|) gives 
the surface state existence criterion 

\ti\<\t2\, (69) 

which coincides with the Shockley criterion ([19]). The 
continuous wave function ([67|) correctly approximates the 
discrete wave function (p!8|) if the decay length is very 
long or equivalently ||ti| — |t2|| <C 1^2 1. However, the 
estimated tunneling amplitudes ti and t2 in Eq. (|6T]) do 
not satisfy the latter condition for Bi2Se3. Therefore, we 
conclude that the discrete Shockley model gives a more 
appropriate description of the surface states in Bi2Se3, 
than a continuous approximation, because the difference 
between \ti\ and \t2\ is rather large. 



1. First-order expansion 



Let us consider the Shockley Hamiltonian (|28|) for the 
fixed value of the in-plane momentum p = p^^ where the 
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FIG. 11: (Color online) Plot of the function t(k) in the com- 
plex plane of (Re t, Im t) . The second-order expansion for t(k) , 
given by Eq. ([77|) with the parameters from Ref. p^] , is plot- 
ted by the solid line for —7v/2 < k < 7r/2. The function t{k) 
in the Shockley model, given by Eq. ()63p with the parameters 
ti and t2 from Eq. (|6ip . is plotted by the dashed line. The 
Shockley contour winds around the origin, which guarantees 
existence of the surface state. 



2. Higher- order expansion 

One may truncate the series for e*^ in Eq. (|64|) at a 
higher order in k 



N 



n=l 



nl 



(70) 



However, such a truncation gives worse continuous de- 
scription of the Shockley surface state. The equation 
t{k) = now has N roots /ci, . . . , /cat. So, there are N 
independent coefficients Cn in a general solution "^{z) = 
Q^^^klz ^ ^ cn^^^^^ to satisfy the boundary condi- 
tions (|65|) and (|66|) . This gives rise to a large number of 
the unphysical surface state solutions, while the Shockley 
model predicts only one surface state. Most of the roots 
kj have large imaginary parts Im kj > 1 . These solutions 
are spurious, because they correspond to the wave func- 
tions decaying over a length shorter than the unit cell of 
the crystal. For example, for TV = 2, Eq. ([7Q|) is 



/2. 



t{k) =ti+t2+it2k 
Then, the equation t{k) = has two roots 

fci,2 = i ± v/-l + 2(l+ti/t2). 



(71) 



(72) 



In the limit |1 + ^1/^2! <C 1, the roots become ki = 
i(l + ^1/^2) and /c2 = 2z. We observe that, while the first 
root ki reproduces the correct approximation Eq. (|68|) . 



the second root k2 has a large imaginary part and must 
be discarded. In another regime, when the expression 
under the square root in Eq. ([72]) is positive, both roots 
have large imaginary parts Im /ci,2 = 1, so the continuous 
approximation is not applicable. Moreover, as pointed 
out in Ref. [34], the continuous description does not dis- 
tinguish between two possible ways of terminating the 
crystal shown in Fig. [TOTa) and (b). A correct bound- 
ary condition should be chosen to distinguish between 
different possible surface terminations. 

Despite these problems, the k'^ terms were kept in the 
effective description of Bi2Se3 in Ref. ^| 



H — Hq + Hi , 

Ho = e{k) ^ {Mo ^ Mie)T, 

Hi = AoTxia X p). 



Bokr, 



OKTy, 



(73) 
(74) 
(75) 



where Mq = -0.28 eV, Mi = 6.86 eVA , Bq = 
2.26 eVA, Ao = 3.33 eVA; Ty and Tx are the Pauli ma- 
trices. In Eq. ([75]) , Hi represents spin-orbit interaction 
and explicitly depends on the spin operators a and the 
in-plane momentum p. Ho depends on the out-of-plane 
momentum k and is responsible for the existence of the 
surface states. Following Ref. [22|, we drop the term e{k) 
in Eq. ([Ti]), because it is proportional to the unit ma- 
trix. Then we apply the unitary transformation e"*'^^^/^, 
which changes —Tx and Tx Tz. So, the Hamilto- 

nian becomes WHoU Ho 



Ho 



Hp) t*{k) \ 
^t{k) -h{p))' 

t{k) = -Mo - Mifc^ + iBok, 
h{p) = Ao{<Txp). 



(76) 

(77) 
(78) 



Now, the Hamiltonian Ho has the same form as in 
Eq. ([28]). Following Ref. we infer that the basis 

for the Hamiltonian ([76]) corresponds to the basis of elec- 
tron orbitals located at the outermost layers of the quin- 
tuplet. Then, the off-diagonal matrix element t{k) in 
Eq. ([77]) corresponds to the second-order expansion of the 
off-diagonal element of the Shockley model ([7T]) , while 
h{p) defines the in-plane dispersion. 

To make explicit correspondence with the previous sec- 
tion, we change units for /c: ka ^ k, where a = 1 nm is 
the size of the elementary cell of Bi2Se3 in the z direction. 
So, we rewrite the parameters Mi/o? Mi, Bo/ a ^ Bo 
in the energy units of eV 

Mo = -0.28 eV, Mi = 0.07 eV, = 0.23 eV. 

Notice that t(/c), Eq. ([77]) , parametrizes a parabola in the 
complex space (Re t,Im t) when k is changed. So we plot 
t{k) defined by Eq. ((TTJ) for -7r/2 < A: < 7r/2 by the solid 
line in Fig.HH Figure [TTI also shows the plot t{k) for the 
discrete Shockley model ([62]) with the parameters ([6T]) by 
the dashed line. We see that the continuous approxima- 
tion to the Hamiltonian agrees with the Shockley model 
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FIG. 12: (Color online) An illustration of the generalized 
Hamiltonian (|8Qp . The Hamiltonian describes a tight-binding 
model with the elementary cell comprised of the A and B 
sites, which are connected via the complex tight-binding am- 
plitudes ti, t2 and ts. 



within a limited range of k with the continuous approxi- 
mation. Nevertheless, the continuous approximation has 
serious deficiencies for construction of the wave functions, 
as described above, whereas the Shockley model provides 
a good overall description for the surface states in Bi2Se3. 



VI. GENERALIZED SHOCKLEY MODEL 

In this section, we generalize the Shockley model to in- 
clude additional inter-cell tunneling amplitudes. To sim- 
plify notations, we present results for the ID case. How- 
ever, the results can be straightforwardly generalized to 
the 3D case by assigning dependence on p = (px^Py) to 
the tunneling amplitudes, as discussed in Sec. [Ill 



A. Additional tight binding amplitude ts 

Let us consider a ID Hamiltonian of the form given by 
Eq. (d]) with 




V 



q 
ts 



(79) 



where the matrix V now contains an additional tight- 
binding amplitude t^. The ID chain model correspond- 
ing to Eq. ([79|) is illustrated in Fig. [121 The amplitude 
ti describes tunneling between the A and B sublattices 
inside the unit cell, and the amplitudes ^2 and ts between 
the unit cells. The introduction of this tight-binding am- 
plitude is motivated by the TI literature [19j, |20|] as well 



by the novel ID models such as the superconducting Ma- 
jorana chain [43, 44] and the Creutz ladder [M^^fiJ- This 
model is a natural mathematical generalization of the 
models considered in the previous sections. The Hamilto- 
nian of the general model has the same form as in Eq. (j5j). 



H{k) 



t*{k) 
t{k) 



with 



t{k) = ti+ t2e'^ + t^e 



(80) 



(81) 



As in Eqs. (p!2]) and (p!3|) , the eigenstate equations for 
the wave functions on the A and B sublattices decouple 
at £^ = 0. (Appendix [Bl proves that the edge state can 
exist only for E = 0.) The zero-energy state on the A 
sublattice has the complex momentum k obtained from 
the equation 



t{k) = ti ^ t2e'^ ^ tse' 



-ik 



0. 



(82) 



We substitute q = e*^ and obtain an equation for the 
rational function t{q) 

t{q) =ti^t2q^tsq-^ =0, (83) 
which has two solutions 



Ql,2 



2t2 ^ 



4t2t3 



(84) 



with the complex momenta ki^2- Using these momenta, 
we construct an edge state that satisfies the boundary 
conditions given by Eqs. (p!Q|) and (pTj) . The edge state 
has the energy Eq = and is localized on the A sublattice 



4-0(2) = 




-En = 0, 



oikiz 



(85) 
(86) 



The wave function (j86]) satisfies the boundary condition 
(pTj) if Im /ci > and Im /c2 > or, equivalent ly. 



\qi\ < 1 and \q2\ < 1. 



(87) 



Likewise, a zero-energy state on the B sublattice has the 
complex momenta k obtained from the equation 



e{k) = tittle 



-ik 



■tie' 



0. 



(88) 



Notice that the symbol of complex conjugation * applies 
only to the tunneling amplitudes in Eq. (|88|) . so t^{k) 7^ 
[t{k)Y if Im/c ^ 0. Equation (j88j) can be obtained by 
replacing k ^ k* in Eq. ([82]) . So, the two solutions /c^ 2 
of Eq. dHHj) and the corresponding q[ 2 can be obtained 
from Eq. ([84]) 



^1,2 ~ ^1,25 ^1,2 ~ 1/^1,2- 



(89) 
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The edge state exists on the B sublattice 



W = ('/i) 




if Im k'i> Q and Im > or, equivalently, 
l^il > 1 and 1^21 > 1. 



(90) 
(91) 

(92) 



To summarize, the edge state (|86]) exists on the A sub- 
lattice if both roots of Eq. ([83|) are inside the unit circle, 
as in Eq. ^ and in Fig. [IS^a). The edge state (|9T]) 
exists on the B sublattice if both roots of Eq. (|83]) are 
outside the unit circle, as in Eq. (|92]) . The edge state 
does not exist if one of the roots is inside and another 
root is outside the unit circle 



l^il > 1 and 1^21 < 1, 



(93) 



as shown in Fig. [TsT b). Obviously, the condi- 
tions ([87j) and (j92j) cannot be met simultaneously, so edge 
state cannot exist on both sublattices simultaneously. 

Like in Sec. Ill A| the criterion for the edge states ex- 
istence can be formulated in terms of the winding num- 
ber of the complex function t{q) along the unit circle 
C = {\q\ = l} 



(94) 



l9l = l 



The criteria given by Eqs. dSZl, and (jMl) are sum- 

marized in the fohowing 



W 



1, edge state ^p^a \z) exists, 

0, edge state does not exist, (95) 

— 1, edge state ^pl°\z) exists. 



To prove it, we use Cauchy's argument principle 
W = Z-P, 



(96) 



which relates the winding number of a complex func- 
tion t{q) on a contour C with the number of zeros Z 
and the number of poles P inside the contour C. Since 
t{q) given by Eq. (|83|) has a pole at g = 0, as shown by 
the thick black dot in Fig. [TSTa) and (b), the number of 
poles is P = 1. The edge state exists on the A sublattice 
if |(7i,2| < 1, in which case W = Z- P = 2-1 = 1. 
The edge state exists on the B sublattice if |gi^2| > 1, in 
which case W = Z — P = — 1 = —1. The edge state 
does not exist for \qi\ > 1 and 1^21 < 1, in which case 
W = Z-P = 1-1 = 0. 

In other words, according to Eq. (j95j), the edge state 
exists if the closed contour 




(b) 
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FIG. 13: (Color online) Topological formulation of the Shock- 
ley criterion. The roots ^1,2 of Eq. (|83p are shown in panels 
(a) and (b) by red dots. An edge state exists if the roots are 
on the same side of the unit circle C = = e*^, /c G (0, 27r)}, 
as shown in panel (a). No edge state exists if the roots are 
on the opposite sides of the unit circle, as shown in panel (b). 
The thick black dot at the origin is the pole of Eq. (|83p . An 
alternative formulation in terms of the winding number (|95p 
is shown in panels (c) and (d). An edge state exists if the 
contour C' = {t{k), k G (0, 27r)} winds around the origin, as 
shown in panel (c). No edge state exists if the contour C' 
does not wind around the origin, as shown in panel (d). 



winds around the origin, as shown in Fig. [TSTc). The di- 
rection of winding of t{k) defines the sublattice on which 
the edge state is localized. An analogous criterion was 
proposed in Ref. [20] (for a comparison with our model, 
see Appendix IC]) . 

When the tunneling amplitudes ti, ^2, and ts are real, 
Eq. (j97j) defines an ellipse 

t{k) = ti + {t2 + ts) cos k + i{t2 - ts) sin k, (98) 

which is shifted by ti from the origin. The ellipse in 
Eq. (|98|) encloses the origin, and the edge state exists if 



1^1 1 < 1^2+^31 



(99) 



Eq. (j99|) represents the generalized Shockley rule of a 
stronger bond: The edge state exists if the broken inter- 
cell bond t2 + ts is stronger than the intra-cell bond ti. 



B. Arbitrary periodic function t(k) 

In the most general form, Hamiltonian (|8Q|) can be 
written as 



C = {t{k), k e (0, 27r)} 



(97) 



H{k) 



r(fe) 

t{k) 



(100) 
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where 

N 

t{k) = ^ tne^'^ (101) 

n=-N 

This model describes a ID tight-binding chain, where 
each unit ceh is coupled to N preceding and N successive 
unit cells. Therefore, for a half-infinite system at z > 1, 
the boundary conditions require that the wave function 
vanishes at the fictitious N sites adjacent to the boundary 

^(-7V + 1) = 0, ^(-1) =0, ^(0) = 0, (102) 

similarly to Eq. ([10]). As in the previous section, let us 
substitute q — e*^ and rewrite Eq. (jlOip in the polyno- 
mial form 

N 

This polynomial has 2N roots. Suppose, the number of 
roots Ni with \qj\ < 1 , j = 1 , • ' ' ? ^ is greater than 
N: Ni > N. In this case, we can construct a trial wave 
function 

^(z) = ^a,^|, (104) 

which vanishes at z oo. The coefficients aj in 
Eq. (|102p are obtained by solving a set of N boundary 
condition equations ()102p . Therefore, in general, there 
are Ni — N linearly independent solutions for the edge 
states localized on the sublattice A. The same result can 
be formulated using the winding number in Eq. ([94]) . In- 
deed, the function (jl03p has a pole of the N-th order at 
q = and A^i zeros at \qj\ < 1, j = 1, . . . , A/"!. Therefore, 
using Cauchy's argument principle (|96|) . we obtain 

W = Z-P = Ni-N. (105) 

Thus, the winding number W of the function t{k) gives 
the number of the edge states. On the other hand, if 
1^ < 0, then there are \W\ degenerate edge states lo- 
calized on the sublattice B. The edge states cannon exist 
simultaneously on both sublattices A and B. There are no 
edge states for 1^ = 0. Finally, in the limit N ^ oo, the 
winding criterion applies to an arbitrary complex func- 
tion t{k) periodic in k. 



VII. CONCLUSIONS 

Like some previous works [l9|, |20|, [32| , our paper ex- 
plores a tight-binding theory of the surface states in topo- 
logical insulators. We show that the surface states can 
be understood using the simple and well-known Shock- 
ley model 113, [in, a ID model with the A and B atoms 
per unit cell, connected via alternating tunneling am- 
plitudes. We generalize the ID Shockley model to the 
3D case described by the 2x2 Hamiltonian (|28|) with 
the diagonal element h{p) and the off-diagonal element 
t{k,p). The diagonal element h{p) deffnes the energy dis- 
persion of the surface states, while the complex- valued 
off-diagonal element t{k,p) defines the domain of exis- 
tence of the surface states. The surface states exist for 
those in-plane momenta p where the phase winding of 
t{k,p) along k G (0, 27r) is non-zero. The sign of the 
winding number gives the sublattice A or B, on which 
the surface states are localized. Equation t{k,p) = 
defines a vortex line in the 3D momentum space 
33], and projection of the vortex line onto the 2D space 
of p is the boundary of the domain where the surface 
states exist. We apply this approach to the TI model on 
the diamond lattice [3]. We show how the evolution of 
the vortex lines is responsible for transitions between the 
"weak" and "strong" TI phases. We discuss why the dis- 
crete Shockley model is better than continuous models 
for the description of the edge states in real materials, 
such as Bi2Se3. The tight-binding model demonstrates 
that different types of surface states are formed depend- 
ing on how crystal is terminated [34]. The surface states 
have the Dirac cone at the center of the Brillouin zone 
when the crystal is cut between the quintuple layers of 
Bi2Se3, but, when the crystal is terminated inside the 
quintuple layer, the surface states have three Dirac cones 
on the boundary of the Brillouin zone. We also general- 
ize the Shockley model to an arbitrary complex function 
t{k) periodic in /c, which includes the long-range inter- 
cell tunneling amplitudes. We prove the validity of the 
winding number criterion in this general case as well. We 
hope that this work will provide a useful toolkit for study- 
ing the surface states in TI, as well as give a transparent 
picture for their physical interpretation. 
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Appendix A: Edge states in the original Shockley 
model 

In Sec. [ni the Schrodinger equation for the wave func- 
tion ^(^) was given in a recursive form for the integer co- 
ordinate z > 1. The main question is whether the recur- 
sion generates a decaying function ^(z) ^ at 2; ^ oo, 
which represents an edge state, or an increasing function 
^(z) ooatz^ 00, which is unphysical. Below, we use 
the generating function method to find convergence cri- 
terion for the edge-state solution "^{z). The Schrodinger 
equation for the original Shockley model ^ is 

V^{z) + (/7 - E)^I/{z + 1) + Vt^(z + 2) = 0, (Al) 

where ^(2:) = [ilja{z)^tljb{z)]^ and ^ > 1, whereas the 
boundary condition is 

{U - E)^{1) + 0^(2) = 0. (A3) 

Let us multiply the z-th Eq. ()Aip by the {z — l)-th power 
of an auxiliary complex variable q and take a sum for 
z > 1 

00 

[V-^{z) + {U - E)-^{z + 1) + + 2)] = 0. 

z=l 

(A4) 

Introducing the generating function 

00 

G(9) = E9"-'*W, (A5) 

Eq. ()A4p can be written as 

[q^V + q{U -E)^ V^]G{q) = V^^{1), (A6) 

where we utilized the boundary condition (|A3p . From 
Eq. ()A6p , we obtain the generating function in terms of 
^(1) 

G{q) = [q^V + q{U - E) ^ V^]'^ ^(1), (A7) 

In order to investigate convergence of ^{z)^ we use the 
following proposition 

Proposition 1. A rational generating function G{q) 
corresponds to an edge state "^{z) - — ^ if and only 
if all poles ^^=1,2,3,... of G{q) have the absolute values 
greater than one, \qj \ > 1. 
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Indeed, a rational function with the poles qj can be 
transformed to the form G{q) = j^i^^^, where fj{q) 
is a polynomial function, and rij is the order of the pole 
qj. Consider a simple example of the first-order pole 
G{q) = = ^z^q/ (liY 1 which corresponds to the geo- 
metric progression. According to Eq. ()A5p . the expansion 
coefficients give the wave function ^f{z) = l/ql~^. Then, 
the absolute values of the pole \qi\ < 1, \qi\ = 1, and 
\qi\ > 1 correspond, respectively, to an un-physical grow- 
ing solution "^{z) ^^^> oo, a bulk state ^(z) = e*^^, and 
a decaying edge state "^{z) - — ^ 0. The case of a more 
complicated G{q) can be reduced to the above simple 
consideration. 

Now let us use Proposition 1 to investigate convergence 
of "^{z). Using Eq. ()A7p and the expressions for U and 
V in Eq. (|M]), we find 



G{q) 



V^a(l) t2 



t2q 



{ti^t2q){t2^tiq)-E^q \ Eq 



. (A8) 



The poles of Eq. (jASp are given by the zeros qi and q2 of 
the denominator, unless they are canceled out by zeros in 
the numerator. Using Vieta's formulas for the quadratic 
equation in the denominator, we obtain ^1^2 = 1- So, if 
qi is greater than one, \qi\ > 1, then q2 is less than one, 
\q2\ < 1. Using Proposition 1, we conclude that there is 
no edge state if the generating function G(q) in Eq. ()A8[) 
has two poles. In order to obtain an edge state, we need 
to reduce the number of poles of the generating function 
G{q). Notice that, if we put E = 0, one pole is canceled 
out, and G{q) greatly simplifies 



G{q) = 



1 



1 + {ti/t2)q 
and the edge state exists if 

h/hl > 1. 



(A9) 



Appendix B: Energy of the edge states in the 
generalized Shockley model 

In this section, we use the generating function method 
to prove that an edge eigenstate for Hamiltonian ([26|) can 
exist only for the eigenenergy E = 0. Like in the previous 
section, Hamiltonian (|26]) can be given in a recursive form 
Eq. (jXTJ) with the following U and V 



V 




U = 



ti 
h 



Using Eq. ()A7p we obtain the generating function 

N{q) 



Giq) = 



D{qy 



(Bl) 



(B2) 



where the numerator 
and denominator 

D{q) = a{q) m - (B4) 
are defined through the polynomials 

m=tlq^^tlq^tl. 
In Eq. (|B3p . the following notation is used for brevity 



(B5) 
(B6) 




(B7) 



According to Proposition 1, the poles of Eq. ()B2p deter- 
mine whether G{q) corresponds to an edge state. The 
potential poles of G{q) are given by zeros of the quar- 
tic polynomial D{q) in the denominator. Thus, let us 
find the structure of zeros of D{q). Suppose, qi is a 
solution of the quartic equation D{qi) = 0. Then, since 
[1^(1/^'*)]* = D{q)/q^ ^ 1/g^ is also a solution of the quar- 
tic equation D{l/ql) = 0. So, in the most general case, 
the polynomial D{q) has zeros qi and ^2, as well as 1/ql 
and 1/^2- Thus, according to Proposition 1, the only 
way to build an edge state is to have the smallest poles 
\qi\ < 1 and \q2\ < 1 canceled out by the zeros of the 
numerator N{q). So, both components of the vector 



N{q) 



tli>2q^ 
tsipiq^ 



■ Expijq 

■ Etp2]q 




(B8) 



must be proportional to {q — qi){q — (72) and, thus, be 
linearly dependent. Hence, the coefficients in front of the 
(AlO) terms q'^ and q^ should also be linearly dependent and so 



0. 



(B9) 



If 1^2! 7^ 1^3 1, then V^iV^2 = 0, so the substitution of 
= and ^^2 7^ (or vice versa) in Eq. (|B8p and 
the requirement, that both components of N{q) are pro- 
portional, lead to ^ = 0. The case 1^2 1 = l^sl is triv- 
ial, because Vieta's formulas for Eq. ()B8p require that 
\qiQ2\ = 1^2/^3! = 1, which contradicts to the initial as- 
sumption that 1^1 1 < 1 and \q2\ < 1. Thus, we have 
proved that the edge states of Hamiltonian ([26|) can only 
exist for £^ = 0, and there are no other edge states. 



17 



Appendix C: Comparison to the model of Ref. [20||. 

Mong and Shivamoggi [l^l considered the tight-binding 
model with the Hamiltonian 

H = Ez [f/*W + V'^iz - 1) + V'l^iz + 1)] , 

(CI) 

Here, U and V represent the intra-cell and inter-cell 2x2 
matrices 

U = rb^, V = rb, (C2) 

where r = {rx^Ty^Tz) are the Pauli matrices in the AB 
sublattice space, bo is a real vector, and b a complex 
vector. Our Hamiltonian ([26]) can also be written in the 



form of Eqs. ([CT]) and ([C2|) with 

b^ = (Reti, Im ti,h), (C3) 
6=(t*+t3,^(t5-t3), 0)/2. (C4) 

Notice that bz = in Eq. ()C4p . It can shown that the 
general Hamiltonian (|C1[) can be always transformed to 
a form with bz = 0. Indeed, the unitary transformations 
^-^rjc/) generated by the Pauli matrices tj rotate the basis 
for the vector b in the bilinear form rb. For an arbitrary 
complex vector 6, it is always possible to select the axis 
z to be orthogonal to both b and 6*, e.g. along i{b x 
6*). So, there always exists a basis where 6^ = 0, and 
the generalized Shockley model discussed in Sec. |VT] is 
equivalent to the model studied in Ref. poj . 



